function [ is ] = SourceMarker( freqlist, magnitude, lN, plotornot )
%SourceMarker plots spectrogram and mark known sources.
%   freqlist  - list of frequencies
%   magnitude - list of magnitudes. mean(magnitudes) ~= 2 !!!
%   plotornot - boolean (0 or 1). Plot spectrogram?
%
%   is        - structure with results
%                 'Name' name of the source
%                 'Frequency' estimated frequency
%                 'Magnitude' estimated magnitude
%                 'Amplitude' estimated amplitude
%                 'HarmonicNumber' harmonic number 
%                                  (root/fundamental harmonic = '1')

    % Convert into colum vectors
    freqlist=freqlist(:);
    magnitude=magnitude(:);
    
    % Known sources
    ks = struct('Name', {}, ...
                'SpinPeriod', {}, ...
                'OrbitalPeriod', {});
            
    % Should be implemented via reading from data file!
    ks(1).objID=1;
    ks(1).Name='Vela X-1';
    ks(1).SpinPeriod=283;
    ks(1).OrbitalPeriod=8.96;
    
    ks(2).objID=2;
    ks(2).Name='GX 301-2';
    ks(2).SpinPeriod=681;
    ks(2).OrbitalPeriod=41.5;
    
    ks(3).objID=3;
    ks(3).Name='A0535+26';
    ks(3).SpinPeriod=104;
    ks(3).OrbitalPeriod=110;
    
    % Plot spectrogram
    if plotornot
        clf
        plot(freqlist, magnitude);
    end
    
    % Local maximums
    lm=zeros(length(freqlist),1);
    lm(2:end-1) = ...
        magnitude(2:end-1)>magnitude(1:end-2) & ...
        magnitude(2:end-1)>magnitude(3:end);
    
    % Local maximum magnitude treshold (>~5 sigma)
    lm = lm & magnitude>12;
    
    % Local maximum frequency treshold (>1e-3 Hz)
    lm = lm & freqlist>1e-3;
    
    % convert lm into list
    lm = find(lm);
    
    % Plot maximums on spectrogram
    hold on
%     plot(freqlist(lm), magnitude(lm), 's');
    
    % Legend list
    legendlist={};
    legendlist{1}='Signal magnitude';
    % Identify
    is = struct('Name', {}, ...
                'objID', {}, ...
                'freq', {}, ...
                'magn', {}, ...
                'hn', {});
            
    for i=1:length(lm)
        % Frequency to investigate
        tempfreq=freqlist(lm(i));
        % Likelihood up to 10th harmonic
        likelihood=zeros(length(ks),10);
        for j=1:length(ks)
            for k=1:10
                temp=tempfreq*ks(j).SpinPeriod/k;
                likelihood(j,k)=max([temp 1/temp])-1;
            end
        end
        % Best candidate
        bestcandidate=find(likelihood==min(likelihood(:)));
        % ID of the source
        sourceid=mod(bestcandidate,length(ks));
        if sourceid==0
            sourceid=length(ks);
        end
        % harmonic number
        harmnum=floor((bestcandidate-1)/length(ks))+1;
        
        is(i).Name=ks(sourceid).Name;
        is(i).objID=ks(sourceid).objID;
        is(i).freq=tempfreq;
        is(i).magn=magnitude(lm(i));
%         is(i).freqer=sqrt(2*lm(i)/lN);
        is(i).hn=harmnum;  
        
        if plotornot
            markerlist='+o*xsd^v';
            plot(freqlist(lm(i)), magnitude(lm(i)), ...
                markerlist(random('unid',8,1,1)),...
                'MarkerEdgeColor', rand(1,3));
            legendlist{i+1}=[is(i).Name ' harmonic #' num2str(is(i).HarmonicNumber)];
        end        
    end
    
    %get rid off dublicates
    marker=ones(1,length(is));
    for i=1:length(is)-1
        for j=i+1:length(is)
            if is(i).objID==is(j).objID && is(i).hn==is(j).hn
                if is(i).magn>is(j).magn
                    marker(j)=0;
                else
                    marker(i)=0;
                end
            end
        end
    end
    is=is(marker==1);
    % Should get rid of ghosts! (leakage)
    
    if plotornot
        legend(legendlist)
    end    
end

